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1  Foreward 

This  report  summarizes  the  research  performed  for  the  Air  Force  Office  of  Scientific  Research 
under  Contract  No.  FA9620-01- 1-0152,  ’’Optimal  Design  of  Smart  Structures.”  The  contract 
began  on  January  1,  2001  and  ran  through  June  30,  2003.  The  Principal  Investigator  was 
Professor  Gordon  G.  Parker,  Department  of  Mechanical  Engineering  -  Engineering  Mechanics 
at  Michigan  Technological  University.  The  co-PI  was  Professor  Bernhard  Bettig,  also  of 
Michigan  Technological  University.  The  AFOSR  program  manager  was  originally  Dr.  Daniel 
Segalman  with  Dr.  Dean  Mook  taking  over  in  2002. 
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2  Objective 

The  objective  of  this  two-year  research  was  to  develop  a  methodology  and  the  required 
analysis  tools  for  simultaneous  optimization  of  active  structures  considering: 

•  the  conventional  structure  topology  . 

•  the  active  material  layout 

•  the  control  system 

This  is  in  contrast  to  the  typical  active  structure  design  approach  of  applying  an  active 
material  treatment  to  an  existing  structure.  Our  hypothesis  was  that  if  the  high  level  design 
goals  were  defined  (e.g.  stiffness,  static  deflection  when  active,  closed-loop  damping  char¬ 
acteristics,  etc.)  without  specifying  the  topology  of  the  structure,  nor  the  active  material 
treatment  and  controller,  then  a  better  final  design  could  be  achieved  by  designing  these 
features  simultaneously  making  the  structure  and  the  controller  work  harmoniously. 
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4  Executive  Summary 

The  primary  objective  of  this  work  was  to  explore  the  hypothesis  that  simultaneous  design  of 
an  active  structure’s  topology,  active  material  treatment,  and  controller  should  yield  better 
performance  than  designing  an  active  material  treatment  for  an  existing  structure,  followed 
by  the  control  design.  Instead  of  tackling  this  problem  in  its  entirety  from  the  onset,  a  step¬ 
wise  approach  was  used.  Specifically,  the  active  structural  design  problem  was  considered 
first  for  static  deflection  and  force  requirement  satisfaction.  This  required  the  development 
of  new  unit  cells,  comprised  of  both  conventional  and  active  material,  for  finite  element 
based  topological  optimization.  After  demonstrating  this  capability  with  several  examples, 
it  was  extended  to  dynamic  analysis.  The  example  considered  was  the  optimal  design  of 
an  active  beam  where  both  the  beam  and  the  active  material  topology  were  optimized 
simultaneously.  This  problem  was  also  decomposed  into  simpler  subproblems  for  initial 
exploration,  specifically 

1.  simultaneous  structure,  sensor,  and  controller  design 

2.  simultaneous  structure,  sensor,  actuator  and  controller  design 

During  this  process  insight  was  gained  not  only  on  cost  function  selection,  but  also  on  how  to 
formulate  the  control  design  portion  of  the  problem.  The  main  conclusions  and  achievements 
were 

•  development  of  a  library  of  finite  element  cells  with  varying  density  of  active  and 
conventional  material  for  use  in  spatial  topological  optimization 

•  development  of  a  virtual  force  approach  for  specifying  structural  strength  requirements 
during  the  optimal  design  process 

•  simultaneous  optimal  design  can  be  used  to  decrease  controller  spillover  effects,  thus, 
increasing  system  performance 

•  topological  optimization  of  both  the  conventional  and  active  material  can  be  used  to 
generate  realizable  designs 

•  the  shape  of  an  active  structure  is  dominated  by  the  shapes  of  the  modes  that  cause 
spillover  instability  when  the  design  goal  is  increased  stability  margin 

Several  unexpected  challenges  arose  during  the  project,  leading  to  additional  interesting 
and  useful  results.  For  example,  to  test  an  optimal  active  beam  design  a  shaped  piezoceramic 
laminate  actuator  was  needed.  Machining  methods  were  investigated  and  it  was  found  that 
a  water  jet  approach  was  suitable  for  cutting,  without  fracturing,  the  ceramic  actuator. 

One  of  the  most  important  aspects  of  this  project  was  the  development  of  an  optimal 
design  code  for  active  structures.  Although  the  development  of  such  a  general  code  was 
not  an  intention  of  the  original  project,  it  seemed  to  be  the  most  logical  approach  as  the 
project  developed.  The  code  uses  a  parallel  (cluster-based)  genetic  algorithm  written  during 
the  project.  Parallelization  was  deemed  necessary  due  to  the  potentially  large  optimization 
problems  that  could  be  encountered  in  the  future.  It  supports  any  user-defined  cost  function 
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and  constraints.  However,  since  it  was  designed  for  structural  optimization,  it  was  linked  to 
ABAQUS  for  general  finite  element  calculations  and  Matlab®  for  rapid  control  design.  Using 
these  off-the-shelf,  general  analysis  codes,  a  wide  range  of  practical  active  structure  problems 
can  be  solved  where  the  system’s  overall  design,  active  material  treatment  and  controller  axe 
designed  together  for  maximum  performance. 
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5  Accomplishments 

In  the  remainder  of  this  report  the  primary  accomplishments  are  documented.  Specifically, 
separate  sections  are  dedicated  to 

•  Development  of  a  new  topological  optimization  unit  cell  that  includes  both  conventional 
and  active  material 

•  An  overview  of  the  optimal  design  code  that  combines  a  parallel  GA  with  ABAQUS 
(for  finite  element  analysis)  and  Matlab®  (for  control  design). 

•  Simultaneous  optimal  design  of  an  active  structure 

5.1  The  Active  Material  Extended  Unit  Cell 

Mechanical  structures  are  designed  to  meet  specific  loading  requirements.  Typically  the  de¬ 
signs  are  iterated  “by  hand” ,  changing  design  features  and  dimensions  again  and  again  until 
a  reasonable  result  is  obtained.  The  question  of  whether  a  chosen  structure  is  optimal  arises 
often.  Topology  optimization,  first  developed  by  Bendspe  and  Kikuchi  (1988)  [12],  automates 
the  process  of  finding  an  optimal  structural  design.  For  a  given  set  of  boundary  conditions 
and  design  specifications  (constraints),  a  topology  is  computed  that  is  optimal  in  terms  of 
a  formulated  cost  function.  A  generic  design  case  is  shown  in  Figure  1  illustrating  several 
boundary  and  loading  conditions.  A  typical  objective  function  is  minimization  of  structural 
compliance  and  a  typical  constraint  is  limiting  the  total  volume  of  structural  material  to  a 
percentage  of  the  design  domain  volume.  Bendspe  and  Kikuchi  defined  the  optimal  struc¬ 
ture  as  the  optimal  distribution  of  material  [12].  To  obtain  different  densities  of  constituent 
material  at  each  point  of  the  design  domain,  a  parametrically  defined  microstructure  with  a 
variable  size  hole  was  introduced.  A  homogenization  technique  was  used  to  obtain  the  effec¬ 
tive  elastic  properties  for  a  material  having  an  infinitely  fine  occurrence  of  the  microstructure 
[3, 18,  4, 10].  To  solve  for  the  optimal  topology  the  design  domain  was  discretized  into  a  grid 
of  finite  elements.  Since  each  finite  element  was  defined  parametrically,  and  could  vary  in 
density  from  0  (no  material)  to  1  (solid  material),  the  topology  optimization  problem  could 
be  transformed  into  a  parameter  optimization  problem. 

This  section  documents  the  development  of  a  new  extended  unit  cell  that  permits  the 
simultaneous  design  of  the  structure  and  its  active  material  treatment.  An  example  is  pro¬ 
vided  for  a  mechanism  that  changes  shape  when  an  electric  field  is  applied  while  achieving 
specified  loading  requirements. 

5.1.1  Homogenization  of  the  extended  unit  cell 

In  order  to  perform  the  topological  optimization,  it  is  necessary  to  know  the  structural 
properties  of  the  unit  cell  as  a  function  of  the  density  parameters  which  axe  selected  during 
optimization.  Using  underlying,  microscopic  structures  to  obtain  these  properties,  as  is 
done  with  homogenization,  results  in  topology  optimization  problems  that  have  proven  to 
be  well-converging  [3,  22],  Compared  to  Simple  Isotropic  Material  with  Penalization  (SIMP) 
approaches  for  obtaining  structural  properties,  homogenization  has  numerical  as  well  as 
physical  advantages. 
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external  forces 


zero  displacement 
in  y*-direction 
no  traction 


Figure  1 :  Generic  design  case  in  topology  optimization. 

Homogenization  is  a  limit  taking  process  as  illustrated  in  Figure  2,  where  a  first  order 
asymptotic  expansion  is  used  to  find  effective  material  properties  [3,  4,  24,  6].  Silva  et  al. 
(1999)  [24]  gives  a  complete  description  of  the  numerical  homogenization  of  single  material 
cells.  A  brief  review  is  given  here  for  the  purpose  of  introducing  the  novel  extended  unit  cell 
in  the  next  section. 
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Figure  2:  Illustration  of  the  homogenization  process. 

The  constitutive  electromechanical  equations  of  a  piezoceramic  material  are: 

&ij  =  Cijkl^kl  ^-kijEki 

A  =  eikieki  +  tikEk- 

where  the  variables  ijkl  take  on  values  1, 2, 3  (1, 2  in  two  dimensional  problems).  The  index 
values  can  also  be  thought  of  as  coordinates  x,  y  and  z  (x  =  1,  y  =  2  and  z  =  3),  where 
Ski  is  the  strain  tensor  and  cry  is  the  stress  tensor.  The  elastic  tensor  is  denoted  by  Cijki- 
Furthermore,  Ek  is  the  electric  field  and  the  electric  displacement  is  denoted  as  A>  while 
denotes  the  dielectric  tensor  and  euj  is  the  piezoelectric  stress  tensor. 

Next,  the  strains  of  equation  (1)  are  expressed  using  the  typical  strain-displacement 
equations 


(2) 
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where  material  displacements  are  denoted  as  it*  in  the  spatial  directions  Xj.  Similarly,  the 
electric  field  is  expressed  in  terms  of  the  electric  field  potential  (f>  as 


d £ 

dx{ 


(3) 


The  differential  form  of  the  mechanical,  dynamical  and  electrostatical  equilibrium  equa¬ 
tions  are  given  by 


da. 


v 


dxj 


+  fi  —0) 


dDj 

dxi 


(4) 


where  fi  are  forces,  and  q  the  charge.  The  constitutive  electromechanical  equations  are 
substituted  into  (4)  yielding  the  weak  form  formulation. 

To  solve  for  homogenized  material  properties,  displacements  axe  expressed  as  the  sum  of 
micro-  and  macro-scale  terms  using  an  asymptotic  expansion  of  the  form  Ui  =  u’-+eu}  +  ... 
and  <j)  =  <f>°  +  e^1  +  . . .  which  are  substituted  into  the  weak  form  formulation.  The  size  of 
the  unit  cell,  denoted  by  e,  is  taken  to  zero  in  the  limit  while  keeping  the  density  constant. 
Following  this,  the  microscopic  problem  is  separated  from  the  macroscopic  problem  yielding 
the  following  equations. 


/(' 


dxlq  dvj  JTWftj 
C£jH  dyi  dyj  +  6klJ  dyk  dyj 


j  [eiU 


dyi  dyi 


tile' 


dT™_  da 
dyk  dy{ 


‘-)  AY  =  / 

jj  J  Vj 

)iY=!e<”i 


d  Y, 


d  Y. 


(5) 


J  ^Cijkl 


dzl  dvi  dZp  dvi 
I"  &kij~ 


dr  =  / 

dm  dVj  '  ~K‘Jdykdyj)  J 

[(  dadzl  dZp  da  \  f 

J  \l  dyt  dyi  eikdykdyi)  J 


e 

'-'pij  i 


ip 


yj 

da 

Vi 


dT. 


(6) 


Mvi,aeH\Y), 

In  these  equations,  p  and  q  take  on  the  same  set  of  values  as  i  and  j  (depending  on  the 
dimensionality  of  the  problem.)  The  domain  of  the  unit  cell  is  denoted  by  Y,  while  and 
a  are  test  functions.  The  functions  xT,  Tpq,  z\  and  2?k  are  characteristic  displacements, 
which  are  calculated  from  the  equations  and  are  constrained  to  be  periodic  in  Y .  The 
calculated  characteristic  displacements  are  then  used  in  the  following  integrations  to  calculate 
the  effective,  homogenized  material  properties. 
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(7) 

(8) 
(9) 

(10) 


To  address  the  active  structure  problem,  a  new  unit  cell  was  developed  that  features  both 
conventional  as  well  as  active  material.  This  extended  unit  cell  is  shown  in  Figure  14(b), 
while  a  classical  unit  cell  featuring  one  material  is  shown  in  Figure  14(a).  The  parameters 
a  and  b  can  vary,  so  that  any  possible  density  of  active  material  (ps),  conventional  material 
(pc)  and  void  (pv)  can  be  obtained,  so  long  as  the  total  volume  of  material  and  void  stays 
constant  (0  <  a  <  b  <  1).  The  densities  of  the  active  and  conventional  materials  are  given 
by  ps  =  1  —  b2  and  pc  =  b2-  a2  respectively.  The  fraction  of  void  per  unit  area  is  denoted  by 
pv  =  a2.  This  enables  the  optimal  material  distribution  task  to  be  formulated  as  a  parameter 
optimization  problem. 


Figure  3:  Unit  cells  (a)  classical  (b)  new  extended  unit  cell. 


To  use  the  homogenization  results  for  topology  optimization  the  homogenization  equa¬ 
tions  of  (7)  are  solved  numerically  for  discrete  values  of  a  and  b.  For  each  combination  of 
values,  the  characteristic  displacements  are  first  calculated  using  a  Galerkin  finite  element 
method.  The  unit  cell  is  discretized  using  quadrilateral,  bilinear  Lagrangian  finite  elements 
with  incompatible  nodes.  The  integrations  are  performed  using  Gaussian  quadrature,  and  a 
sparse-matrix  based  inversion  method  is  used  to  solve  the  resulting  set  of  linear  equations. 
The  finite  element  formulation  features  two  mechanical  degrees  of  freedom,  and  one  electrical 
degree  of  freedom  per  node.  There  are  four  nodes  per  finite  element.  Plane  strain  elements 
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are  used,  corresponding  to  a  very  thick  plate.  Further  details  regarding  the  implementation 
of  the  numerical  method  are  available  in  Buehler  et  al.  [6]  and  Bathe  [2].  These  discrete  ho¬ 
mogenization  results  are  then  used  to  perform  a  polynomial  surface  fit  yielding  an  analytical 
function  for  each  material  coefficient.  A  6th  order  polynomial  was  found  to  be  sufficient. 

Figures  4  and  5  show  typical  results.  In  this  example,  the  piezoelectric  material,  by  itself, 
has  elasticity  coefficients  am  —  2.1978  x  107,  C1122  =  0.6593  x  107  and  ci 212  =  0.7692  x  107. 
The  piezoelectric  stress  coefficients  are  e2 11  =  6122  =  0.1  and  the  permittivity  is  = 
1.140  x  10-7.  The  conventional  material,  by  itself,  has  elasticity  coefficients  am  =  1.12  x  107, 
C1122  =  0.323  x  107  and  cm2  =  0.399  x  107.  Both  materials  axe  isotropic.  The  homogenized 
material  elastic  properties  Cijki  are  plotted  in  Figure  4.  The  effective  piezoelectric  stress 
coefficient  e\22  =  e211  and  permittivities  =  e%2  as  a  function  of  the  parameters  a2  and  b2 
are  shown  in  Figure  5,  where  all  quantities  are  assumed  to  be  unit-less.  The  behavior  of  the 
composite  approaches  the  pure  material  when  a  =1  and  b  =  0,  or  vice  versa.  All  material 
coefficients  not  shown  are  zero. 


FIGURE  4:  The  effective  elastic  coefficients  cf>-kl  as  a  function  of  a2  and  b2. 

The  piezoelectric  material  we  use  is  a  “fictitious”  isotropic  piezoelectric  material  without  a 
poling  direction.  In  “classical”,  “real”  active  materials,  there  usually  exists  a  poling  direction, 
for  example  that  determines  in  which  orientation  the  material  expands  when  an  electrical 
field  is  applied.  The  concept  of  using  this  material  model  is  motivated  by  the  desire  to  have  a 
material  that  expands  in  the  y  direction  once  an  electrical  field  is  applied  in  the  x  direction, 
and  expands  in  the  x  direction  once  an  electrical  field  is  applied  in  the  y  direction.  The 
expansion  is  always  orthogonal  to  the  electric  field  vector.  Once  the  electric  field  vector  is 
known,  the  model  piezoelectric  material  can  be  replaced  by  real  piezoelectric  material  which 
has  the  poling  direction  according  to  the  electric  field  vector. 
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FIGURE  5:  The  effective  piezoelectric  stress  coefficient  ef22  —  e2n  an<^  effective  permittivity  coef¬ 
ficients  ejj  =  e22  as  a  function  of  a 2  and  b2. 


In  the  discretization  of  the  design  domain  into  finite  elements,  a  different  microstructure 
is  admitted  in  each  element  leading  to  different  densities  in  each  element.  The  whole  set 
of  design  parameters  is  denoted  by  a  =  (a,  b).  The  expression  of  material  properties  as  a 
function  of  (a,  b)  makes  it  possible  to  treat  the  microstructure  by  effective  material  properties, 
instead  of  modeling  the  microstructure  completely. 

5.1.2  Design  of  a  Structure  and  Embedded  Actuator  with  Maximum  Activated 
Deflection 

The  goal  is  to  synthesize  a  structure  that  converts  electrical  into  mechanical  energy  where 
the  cost  function  is  defined  by  the  displacement  of  a  specific  node.  With  this  scheme,  smart 
structures  are  designed  that  output  a  deflection  at  a  defined  location  when  electrical  power 
is  applied.  In  this  example,  the  tip  deflection  of  the  lower  left  node,  denoted  (xo,yo),  as 
shown  in  Figure  6  is  to  be  maximized.  We  assume  that  a  constant  electric  field  is  applied  in 
the  design  domain  with  a  constant  gradient  of  the  electric  potential  in  the  direction  shown. 
In  addition  to  the  pure  displacement-based  cost  function,  we  introduce  a  stiffness  criterion 
using  an  additional  virtual  boundary  force  Fv  acting  in  the  opposite  direction  of  the  desired 
displacement.  The  criterion  does  not  appear  in  the  cost  function  or  the  constraints  directly, 
but  is  only  present  as  a  boundary  condition  when  the  finite  element  problem  is  solved. 

Since  the  electric  field  is  only  in  one  direction,  the  “fictitious”  material  is  only  actuated 
in  one  direction,  making  it  equivalent  to  an  e31  material  (actuated  in  the  up  and  down 
direction).  At  this  point,  the  electric  field  is  just  given  without  any  consideration  for  how  it 
can  be  realized.  This  is  discussed  later  in  the  discussion  section. 

Usually,  in  compliant  mechanism  design,  or  piezoelectric  actuator  design,  optimization 
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FIGURE  6:  Design  case  for  maximizing  the  tip  deflection  with  the  virtual  force  for  the  structural 
stiffness  criterion. 


of  displacement  of  a  given  point  is  achieved  by  maximizing  the  mean  mutual  energy,  or  the 
mean  transduction  [19,  25,  21].  To  provide  stiffness,  the  mean  compliance  is  minimized  in 
both  cases.  In  this  example,  we  propose  a  somewhat  different  scheme  to  achieve  the  same 
goal,  which  is  large  displacement  on  the  one  hand,  and  sufficient  stiffness  of  the  structure 
on  the  other  hand.  Our  approach  of  optimizing  the  deflection  at  a  chosen  point  is  related 
to  the  maximization  of  an  equivalent  mutual  energy  (or  mean  transduction).  However,  the 
stiffness  criterion  is  achieved  by  applying  a  dummy  load. 

The  stiffness  criterion  and  the  deflection  criterion  are  essentially  counteracting,  and  a 
tradeoff  is  required  in  the  final  solution.  In  very  compliant  structures,  one  can  have  large 
displacements,  but  no  force.  On  the  other  hand,  in  very  stiff  structures,  there  is  little 
displacement,  but  large  forces.  In  the  classical  approaches,  prior  to  optimization  a  parameter 
has  to  be  chosen  which  determines  the  relative  strength  of  the  stiffness  criterion  versus  the 
deflection  criterion  [19,  25,  21].  In  the  method  illustrated  in  this  example,  the  relative 
strength  has  to  be  chosen  carefully  to  obtain  useful  designs,  here  by  selecting  the  strength 
of  the  force  versus  the  strength  of  the  electric  field. 

The  cost  function  is  J( a)  =  -ux(xo,yo),  capturing  the  x  displacement  of  the  lower  left 
node.  The  displacements  ux  are  subject  to  the  “constraints”  of  satisfying  the  equilibrium 
equations  of  (4).  These  are  implemented  by  defining  four  symmetric  bilinear  operators  (for 
details,  see  e.g.  [16]): 


,  ,,  f  dukdvi 

i(u,v)  =  J 

f  dm  da 

c(u,a)=  J  e^^-dn, 


b{(j>,v)  =  / 

‘  9(f)  dvk 

&ikl  &  o  Cli  Z, 
OXiOXi 

d(<f>,  a)  =  J 

f  d<f>  da  JO 

€ik  o  rv  df2, 

UX%  (sXk 

Q 


where  Vi  and  a  are  test  functions  that  satisfy  the  boundary  conditions.  Note  that  two 
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operators  are  equivalent  c(u,  a)  =  b(a,u).  The  right  hand  sides  of  the  piezoelectric  equations 
of  (4)  are  represented  by  four  additional  operators, 

Fm;r(t,  v)  =  J  t&i  dT,  (12) 

r 

Fe;r(q,a)  =  J  qa  dr,  (13) 

r 

Fmln(F,v)  =  /  FiVidQ,  (14) 

si 

Fe;o(Q,a)  =  J  QadO,  (15) 

SI 


where  L  are  the  traction  forces,  q  are  the  applied  charges  at  the  boundaries,  F  and  Q  are 
the  volumetric  forces  and  charges  respectively.  This  allows  the  piezoelectric  equations  to  be 
expressed  compactly  as 

a(u,  v )  +  b(<p,  v )  =  Fm;r(t,  v)  +  F  m-,si(F,  v), 
c(u,  a)  —  d(<f>,  a)  =  Fe;r(g,  a)  +  Fe;n(Q,  a), 


Vwj.aGfT1^), 

where  O  refers  to  the  design  domain.  The  optimization  problem  can  then  be  summarized  as 
finding  a  set  of  elemental  design  parameters,  a  such  that 


min  :  J(a)  —  —ux, 

a>i,bi 

subject  to  :  a(u,  v )  +  b(<p,  v )  =  F v )  +  Fm;n (F,  v), 
c(u,  a)  -  d((f>,  a)  =  Fe;r (Q,  a)  +  Fe;f2(g,  a), 

Vvi,  a  e  ifx(0),  J  pc dft  =  VCV0,  J  PsdCt  =  VSVQ, 

o  h 

0  <a,i<  1,  di  <bi  <  1. 


(17) 


It  should  be  noted  that  constraints  are  also  imposed  on  the  amount  of  material  and  on  the 
elemental  design  parameters  as  described  by  equation  (17).  Specifically,  volume  constraints 
on  active  and  conventional  material  are  imposed,  denoted  by  Vc  and  ,  where  Vo  is  the  area 
of  the  design  domain.  In  addition  the  design  parameters  are  required  to  be  continuous  and 
between  zero  and  unity. 
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5.1.3  Solution  of  Optimization  Problem 

Sequential  Quadratic  Programming  (SQP)  is  used  to  solve  the  optimization  problem  of 
equation  (17).  At  each  iteration  the  cost  function  J(A)  is  obtained  using  the  finite  element 
method.  In  the  piezoelectric  case,  in  addition  to  the  stiffness  matrix  Kuu,  there  are  also  the 
piezoelectric  matrix  Ku^,  as  well  as  the  dielectric  matrix  K^. 

The  displacement  of  all  nodes  is  in  general  given  by  u  =  K-1F,  where  K  is  the  assembled 
stiffness  matrix,  F  is  the  force  vector,  and  u  is  the  displacement  vector.  It  should  be  noted 
that  both  K  and  F  depend  on  the  design  parameters  a.  A  constant  electrical  field  is  applied 
by  specifying  the  electrical  potential  throughout  the  design  domain.  The  elemental  force 
vector  is  therefore  given  by  F  =  -Ku^  .  In  the  present  case  for  a  constant  electrical  field 
in  the  y  direction,  the  vector  (f)  =  [-</>0, 0, 0,  -<f> 0]T  (4>o  is  related  to  the  strength  of  the 
electric  field,  according  to  equation  (3)).  Therefore,  the  problem  can  essentially  be  solved  as 
a  conventional,  completely  elastic  problem. 

The  cost  function  of  (17)  can  be  expressed  using  an  additional  mapping  vector  B,  that  is, 
J{a)  =  B(r)K-1(a)F(a).  With  this  formulation  any  one  displacement  can  be  used  as  the 
cost  function.  Here,  we  simply  need  B  =  [0  ...  0  1  0].  Since  the  SQP  method  is  gradient 
based,  we  also  need  the  partial  derivative  of  J(a)  with  respect  to  the  design  parameters. 
This  can  be  obtained  directly,  in  the  general  form,  from  the  finite  element  analysis  results  as 

=  -B(r)K"1(a)^|^K-1(a)F(a)  +  B(r)K"1(a)^|^.  (18) 

We  note  that  since  the  optimization  problem  is  nonlinear,  local  minima  can  occur.  A 
question  is  therefore  whether  the  obtained  solution  is  a  local  or  a  global  minimum.  To 
address  this  the  optimization  problem  was  solved  for  several  different  initial  a  distributions. 
In  each  case,  similar  results  were  obtained  providing  confidence  that  the  optimal  distribution 
of  smart  and  conventional  material  is  close  to  a  global  minimum  of  the  constrained  objective 
function,  however  this  is  not  intended  to  be  proof  of  global  optimality. 

5.1.4  Results 

In  the  first  set  of  results  a  finite  element  discretization  of  40  x  20  was  used  in  the  topology 
optimization  problems.  The  “raw” ,  unprocessed  optimization  results  for  two  discretizations 
in  Figure  7(a)  for  an  active  material  volume  constraint  of  20  percent,  and  a  volume  constraint 
of  30  percent  for  conventional  material.  In  Figure  7(f)  for  an  active  material  volume  con¬ 
straint  of  30  percent,  and  a  volume  constraint  of  40  percent  for  conventional  material.  Little 
checkerboarding  is  observed  in  both  examples,  and  the  smart  and  conventional  material  is 
well-separated. 

It  is  often  advantageous  to  post-process  the  data  further.  This  is  illustrated  in  Figure  7(b) 
and  Figure  7(g).  A  finite  element  was  considered  comprised  of  completely  active  material 
if  ps  >  0.50.  Similarly,  an  element  was  considered  comprised  of  completely  conventional 
material  when  pc  >  0.50.  Finite  element  analysis  results  are  shown  in  Figure  7(c)  and  (h). 
An  applied  force  leads  to  a  deformation  of  the  structure,  as  shown  in  Figure  7(c)  and  (h)  with 
no  field  applied.  The  electric  field  is  turned  on  in  Figure  7(d)  and  (i)  illustrating  the  negation 
of  the  force  induced  displacement.  Figure  7(e)  and  (j)  shows  a  possible  interpretation  of  the 
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FIGURE  7:  Optimized  designs  for  two  sets  of  volume  constraints.  Conventional  material  is  shown 
black  and  active  material  is  shown  grey.  From  left  to  right:  (a,f)  raw  optimization 
results,  (b,g)  interpretation  of  the  results,  (c,h)  deflection  with  force  applied  and  elec¬ 
trical  power  switched  off,  (d,i)  negation  of  deflection  when  electrical  power  is  switched 
on,  (ej)  possible  “continuous”  interpretation. 

resulting  design  using  a  contiguous  segment  of  active  material.  An  interesting  observation  is 
that  it  seems  to  put  some  of  the  conventional  material  in  a  place  where  it  doesn’t  do  anything 
(observe  the  vertical  strip  on  the  right  in  Figure  7  (a)-(d)).  If  this  material  were  used  on  the 
left,  likely  it  would  make  the  structure  stiffer,  resulting  in  a  smaller  displacement  when  the 
structure  is  activated. 

Another  set  of  results  of  topology  optimization  of  a  clamp  mechanism,  systematically 
changing  the  volume  constraint  for  the  active  material  and  the  conventional  material,  as 
well  as  changing  the  shape  of  the  design  domain  from  rectangular  (1)  to  square  (2)  is  shown 
in  Figure  8.  In  Figures  8(a)-(c),  the  active  material  constraint  is  kept  constant  at  30  percent, 
while  the  conventional  material  volume  constraint  varies  from  8  to  30  percent.  In  Figures  8(c) 
and  (e),  the  active  material  constraint  is  varied  from  30  percent  to  8  percent,  while  the  volume 
constraint  for  active  material  is  kept  constant  at  30  percent.  The  smart  material  tends  to 
be  put  close  to  the  fixed  boundary.  The  optimization  runs  also  reveal  that  the  shape  of  the 
conventional  material  changes,  when  the  amount  of  smart  material  is  modified. 

The  resulting  design  is  obviously  affected  directly  by  the  volume  constraints,  however, 
the  magnitude  of  the  virtual  force  acting  on  the  displaced  node  must  also  be  chosen  carefully. 
By  changing  the  strength  of  the  force,  it  is  possible  to  affect  the  stiffness  of  the  structure. 
Figure  9  (a)  shows  the  resulting  structures  for  different  strengths  of  the  virtual  force.  In 
this  figure  the  gray  colored  shapes  (deflecting  to  the  left)  show  the  structures  with  just  the 
deflection  due  to  the  piezoelectric  effect.  The  darker  (colored)  shapes  show  the  structures 
with  the  piezoelectric  effect  and  with  the  application  of  a  constant  force.  In  panels  1-4,  the 
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Figure  8:  Topology  optimization  of  a  clamp  mechanism:  systematically  changing  the  volume 
constraints,  as  well  as  changing  the  shape  of  the  design  domain  from  rectangular  (1)  to 
square  (2). 

virtual  force  is  decreasing.  It  can  be  observed  that  if  the  force  is  small  the  displacement  will 
be  large  when  the  PZT  is  activated,  but  the  structure  will  be  very  compliant  and  will  not 
be  able  to  support  even  the  small  virtual  force.  It  is  also  observed  that  if  the  virtual  force  is 
small,  the  structures  often  feature  hinges.  On  the  other  hand,  if  the  virtual  force  is  large  the 
structure  will  be  stiff  and  the  resulting  displacements  will  be  small,  both  with  and  without 
the  virtual  force.  Figure  9  (b)  shows  the  tip  deflection  with  and  without  a  force  applied,  as 
a  function  of  the  virtual  force  applied  during  the  optimization. 

The  results  in  Figure  9  also  indicate  that  there  exists  a  convergence  to  an  optimal  solution 
once  the  virtual  force  is  below  or  above  a  critical  value  (below  6,000  the  resulting  tip  deflection 
does  not  change  much,  as  for  magnitudes  of  the  virtual  force  above  11,000).  Such  behavior 
may  be  found  in  many  design  cases.  However,  there  is  no  indication  that  the  actual  values 
would  be  identical  in  different  problems.  Therefore,  it  seems  necessary  that  different  values 
should  be  tried  for  a  specific  problem  and  then  decided  what  is  most  suitable  to  achieve  good 
results. 

Finally,  we  discuss  the  convergence  of  the  different  problems  studied  in  this  work.  Fig¬ 
ure  10  (a)  plots  the  development  of  the  objective  function  over  76  iterations,  for  the  four 
different  optimization  cases  shown  in  Figure  9.  It  can  be  seen  that  the  optimum  solution  is 
approached  monotonically  in  each  case. 

In  our  calculations,  the  optimization  is  stopped  when  the  material  distribution  is  such 
that  densities  close  to  zero  or  one  dominate  in  the  design  domain.  This  criterion  could 
further  be  used  to  define  a  critical  slope  of  the  objective  function  with  respect  to  iterations 
A  J.  The  optimization  is  then  stopped  once  the  change  in  cost  becomes  sufficiently  small, 
e.g.  below  a  threshold  of  A  J  =  0.01  between  two  iteration  steps. 
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FIGURE  9:  Varying  the  strength  of  the  virtual  force.  Gray  structure  shows  the  results  when  only 
the  electric  field  is  applied.  Dark  (color)  structure  shows  the  results  when  a  force  of 
constant  magnitude  is  applied.  Graph  shows  the  tip  deflection  with  and  without  the 
force  applied,  as  a  function  of  the  virtual  force. 


It  is  interesting  to  note  that  when  a  larger  virtual  force  is  used,  the  convergence  is  also 
faster  (cf.  cases  1  and  2  with  large  virtual  force,  and  cases  3  and  4  with  smaller  virtual 
force). 

Figure  10  (b)  plots  the  constraint  errors  over  the  iterations  for  the  four  optimization 
cases  shown  in  Figure  9.  Figure  10  (b)  proves  that  the  constraints  are  always  satisfied 
within  an  average  deviation  below  1  x  10-4,  at  least  for  the  cases  1-3  (when  the  virtual  force 
gets  very  small  as  in  case  4,  convergence  becomes  poor).  Convergence  in  the  other  cases 
considered  (Figures  7  and  8)  is  very  similar  to  the  results  shown  in  Figure  10  (a),  case  1  and 
2.  Overall  the  optimizations  typically  required  60-80  iterations  to  converge  while  satisfying 
the  constraints  below  an  error  of  10-4. 

5.1.5  Discussion 

The  computations  show  that  it  is  feasible  to  perform  the  numerical  homogenization  of  a  unit 
cell  with  parametrically  defined  densities  of  active  material,  conventional  material,  and  void. 
They  also  show  that  it  is  feasible  to  use  the  resulting  homogenized  properties  to  perform  a 
topological  optimization  in  which  there  are  no  constraints  on  the  position  of  the  conventional 
or  active  material.  As  well,  it  was  found  that  the  different  constituting  materials  are  well- 
separated  in  the  final  solution. 
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Figure  10:  Development  of  the  objective  function  (a)  and  the  constraint  error  (b)  over  the  itera¬ 
tions,  for  the  three  different  optimization  cases  shown  in  Figure  9. 


We  find  that  it  is  crucial  to  include  a  stiffness  criterion  in  order  to  obtain  good  designs, 
as  clearly  shown  in  Figure  9.  This  is  in  agreement  with  previous  research  on  the  topic 
[19,  25,  21].  Further  discussion  on  the  topic  of  competing  design  for  compliance  and  stiffness 
for  conventional  structures  could  be  found  in  a  paper  by  Chen  et  al.  (2001)  [8]. 

It  was  observed  that  for  very  low  volume  fractions  of  both  acvtive  and  conventional 
material,  the  method  is  not  as  effective  as  in  other  cases.  However,  as  can  be  seen  in 
Figure  7(c),  even  for  low  volume  fractions  of  active  material  good  convergence  is  obtained. 

The  material  we  used  in  the  examples  is  characterized  by  its  piezoelectric  coupling  co¬ 
efficients  e2n  and  e122,  while  all  other  piezoelectric  stress  coefficients  were  assumed  to  be 
zero.  These  piezoelectric  properties  were  chosen  in  order  to  achieve  a  piezoelectric  effect 
orthogonal  to  the  applied  electric  field.  Other  choices  (e.g.  em  and  6222)  axe  possible  as 
well. 

The  principle  aim  of  this  research  was  not  to  synthesize  structures  that  axe  manufac¬ 
turable,  but  to  develop  general  and  basic  principles.  In  terms  of  manufacturing  the  resulting 
structures,  the  proposed  material  microstructure  is  not  literally  acceptable.  But,  as  the  re¬ 
sults  show,  regions  with  either  one  or  the  other  material,  or  void,  appear  in  the  solution.  As 
a  result,  manufacturing  of  a  mixture  of  materials  is  not  required. 
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The  resulting  topologies  are  also  post-processed  in  order  to  interpret  realizable  structures 
that  can  be  built.  An  interpretation  using  stack  actuators  is  proposed  which  are  added  to 
the  conventional  structure  in  such  a  way  that  they  have  the  desired  effect.  The  electric 
field  leads  to  an  expansion  (respectively  contraction)  in  the  orthogonal  direction.  Analyzing 
the  resulting  displacement,  one  can  orient  stack  actuators  (e.g.  3-3  piezoelectric  material) 
properly  and  overcome  the  difficulty  of  having  a  2-1  material  with  planar  electric  fields. 
Generally,  if  the  active  material  covers  the  entire  cross-section  of  a  slender  member  such  as 
the  red  area  in  Figure  7(e)  and  (j),  then  it  it  may  be  possible  to  replace  that  with  stack 
actuators.  However,  if  only  a  part  of  the  member  cross  section  is  made  of  active  material, 
then  under  electric  forces  the  member  would  be  bending  and  it  may  not  be  possible  to 
directly  come  up  with  a  manufacturable  part  that  will  replace  that  member. 

In  some  cases  this  allows  building  the  resulting  structures.  For  instance,  the  part  in 
Figure  7(e)  contains  active  material  (red  region)  that  could  be  realized  by  many  stack  layers. 
Using  different  sizes  of  these  piezoelectric  stack  actuators,  oriented  in  the  proper  orientation, 
even  the  optimal  shape  could  be  approximated.  In  this  case,  the  electrodes  would  be  oriented 
orthogonal  to  the  direction  of  the  active  material  beam  portion  of  the  topology  optimization 
result. 


5.2  Parallel  Genetic  Algorithm  Active  Structure  Design  Code 

Both  sequential  quadratic  programming  (SQP)  and  genetic  algorithm  (GA)  methods  were 
used  to  solve  the  optimization  problems  in  this  project.  Specifically,  the  fmincon  function 
in  Matlab®  ,  an  SQP  code,  was  used  to  solve  low  dimensional  optimization  problems  con¬ 
taining  no  more  than  50  design  parameters.  This  was  the  case  for  the  static  design  example 
considered  in  Section  5.1  and  for  the  optimal  sensor  design  study  that  is  described  in  Section 
5.3.1.  Although  these  problems  were  relatively  small,  it  was  found  that  running  Matlab®  on 
a  single  computer  would  be  problematic  for  larger  optimization  problems  and  those  requiring 
complex  cost  function  evaluations. 

Based  on  this  realization  it  was  decided  that  a  general,  parallel  processing  based  approach 
should  be  developed.  A  genetic  algorithim  method  was  eventually  selected  since  GAs  can 
accomodate  both  discrete  and  continuous  cost  functions.  A  computing  cluster  formulation 
was  chosen  since  these  are  more  readily  available  to  researchers  and  are  scalable,  as  compared 
to  single,  multiprocessor  computers. 

GA  Code  Description  The  GA  code  was  written  in  C  using  the  Message  Passing  Inter¬ 
face  (MPI)  standard  for  sharing  data  between  computers  in  the  cluster.  At  a  high-level,  a 
GA  code  starts  with  an  initial  population  consisting  of  N  members.  Members  are  defined 
by  a  genetic  code  that  maps  to  discrete  values  of  the  design  parameters.  After  initializing 
the  population  with  randomly  created  members,  all  of  their  cost  functions  are  computed  and 
stored.  New  generations  are  created  by  cross-over  and  mutation,  where  the  cost  function  of 
each  new  member  must  again  be  computed.  When  a  new  generation  is  created,  the  lowest 
cost  function  members  are  removed.  This  process  continues  for  a  fixed  number  of  genera¬ 
tions,  or  until  some  stopping  criteria,  for  example  based  on  the  change  of  the  member’s  net 
cost,  is  satisfied.  A  flow  chart  of  the  GA  code  developed  during  this  project  is  shown  in 
Figure  11. 
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At  a  lower-level,  the  parallel  implementation  allowed  the  available  computers  in  the 
cluster  to  work  on  computing  the  cost  functions  for  the  members,  thus  reducing  the  time 
needed  to  evaluate  a  new  generation.  A  ‘master’  computer  was  given  the  task  of  organizing 
how  the  members  were  distributed  amongst  the  cluster.  The  GA  code  was  written  with 
sufficient  generality  to  accomodate  any  number  of  computers  in  the  cluster  and  was  tested 
with  up  to  47  Sun  workstations. 

One  important  benefit  of  the  GA  approach  over  gradient  based  methods  is  that  a  set  of 
near-optimal  solutions  is  provided  that  may  be  in  vastly  different  places  in  the  design  space. 
This  provides  the  designer  a  set  of  different  designs  to  choose  from,  facilitating  the  crucial 
capability  of  bringing  human  decision  making  into  the  optimization  process.  To  implement 
this  in  a  practical  way,  the  capability  was  added  to  the  code  allowing  a  user  to  look  at 
the  best  designs  as  the  code  is  running.  Additionally,  the  user  can  modify  aspects  of  those 
designs,  and  re-insert  them  into  the  population.  This  real-time  ’’seeding”  capability  is  a 
unique  feature  of  the  design  environment  and  is  intended  to  bring  together  some  of  the  best 
aspects  of  human-based  design  and  computer-based  design. 

Multidiscplinary  Design  Environment  There  are  three  primary  tasks  for  simultaneous 
active  structure  design 

1.  selection  of  parameter  values  that  uniquely  define  the  system  consisting  of  both  struc¬ 
tural  and  control  attributes 

2.  creation  and  analysis  of  the  corresponding  static  or  dynamic  model  defined  by  the 
structural  parameters  selected  in  step  1. 

3.  development  of  the  control  law  based  on  the  model  created  in  step  2  and  the  control 
law  parameters  selected  in  step  1. 

In  the  most  general  case  it  was  expected  that  finite  element  analysis  would  be  required 
so  that  mode  shape  and  frequency  information,  a  model  of  the  structure,  would  be  available 
for  the  control  law  design  phase.  For  this  purpose  ABAQUS,  a  commercially  available  fi¬ 
nite  element  code,  was  chosen.  Similarly,  to  encompass  the  most  general  control  law  design 
possible,  it  was  decided  to  facilitate  the  automatic  inclusion  of  Matlab®  into  the  optimal 
design  code.  Using  low-level  communication  protocols  within  both  ABAQUS  and  Matlab® 
it  was  eventually  possible  to  create  a  fully  integrated  system  where  the  GA  generated  the 
design  parameters,  passed  the  structural  information  to  ABAQUS,  along  with  the  analysis 
instructions.  The  ABAQUS  results  were  then  automatically  passed  to  Matlab®  ,  along  with 
control  design  instructions  and  the  control  related  parameters  from  the  GA.  With  the  cost 
function  evaluation  definition  orginating  from  the  GA,  Matlab®  could  then  design  the  con¬ 
troller  and  generate  the  cost  function  based  on  either  time  domain  simulation  or  frequency 
domain  analysis.  A  block  diagram  illustrating  the  process  is  shown  in  Figure  12. 

This  general  framework  was  certainly  sufficient  for  the  active  material  systems  considered 
in  the  examples  of  the  following  section.  However,  it  has  a  broader  range  of  applicability 
including  automatic  design  of  active  truss  systems,  wing  internal  structures,  etc.  Other 
application  areas  are  currently  begin  explored  with  Air  Force  agencies  and  DARPA. 
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FIGURE  1 1 :  Flow  of  the  genetic  algorithm  based  active  structure  design  code.  Green  blocks  indicate 
where  fitness  function  calculations  are  required  and  are  shown  in  detail  in  Figure  12 
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5.3  Simultaneous  Optimal  Design  of  an  Active  Structure 

As  described  in  the  Executive  Summary,  an  incremental  approach  was  used  leading  up  to 
full  simultaneous  design  of  a  structure’s  conventional  and  active  material  as  well  as  the 
controller.  The  results  of  these  two  steps  are  documented  in  the  remainder  of  this  section. 
First  the  simultaneous  design  of  the  sensor  and  controller  package,  given  an  existing  structure 
is  described.  This  exercise  provided  valuable  insight  into  the  cost  function  selection,  and  the 
approach  for  control  design.  Second  the  full  simultaneous  design  problem  is  considered  in 
Section  5.3.2. 

5.3.1  Simultaneous  Optimization-of  the  Sensors  and  Controller 

Simultaneous  optimization  of  the  smart  material  distribution  and  the  control  strategy  was 
investigated  with  the  goal  of  improving  the  stability  margin  of  a  closed-loop  active  struc¬ 
ture.  This  required  extending  the  topological  optimization  for  static  structure  requirements, 
described  in  section  5.1,  to  the  dynamic  case.  In  addition,  investigation  of  suitable  cost 
functions  was  performed.  Through  topological  optimization  of  the  active  material  sensor, 
using  a  homogenization  approach  and  a  linear  quadratic  regulator  (LQR),  a  new  type  of 
sensor,  with  the  ability  to  increase  the  closed-loop  stability  margin,  was  obtained. 

One  constraint  that  must  be  specified  at  the  beginning  of  any  topological  optimization 
problem  is  the  design  domain  size.  In  section  5.1  this  was  addressed  in  terms  of  volume 
fraction  of  active  versus  conventional  material.  In  this  section  two  different  sensor  design 
domains  are  considered  (1)  discrete  (2)  segmented  distributed  as  shown  in  Figure  13.  This 
type  of  comparison  was  performed  to  provide  some  insight  into  how  the  sensor  design  domain 
should  be  specified  for  a  general  structure.  The  primary  conclusion  was  that  the  segmented 
distributed  domain  provided  better  closed-loop  performance,  from  a  stability  margin  perspec¬ 
tive. 

Introduction  A  controller  that  satisfies  the  design  requirements  for  all  admissible  pertur¬ 
bations  is  termed  robust[7].  As  applied  to  the  control  of  flexible  structures,  this  means  that 
a  closed-loop  controller  must  be  able  to  operate  in  the  presence  of  modelling  uncertainties 
that  lead  to  spillover  instabilities.  This  has  been  studied  extensively  [1][15][17]. 

In  this  approach,  a  method  for  sensor  design  that  increased  the  stability  margin  by  simul¬ 
taneously  optimizing  a  Polyvinylidene  Fluoride  (PVDF)  sensor  distribution  and  the  control 
strategy  for  a  closed-loop  system  was  developed.  The  approach  was  motivated  by  Collet  [9] 
who  illustrated  optimal  sensor  design  based  on  minimizing  the  observability  gramian  of  the 
open  loop  system.  Here  a  more  general  sensor  design  strategy  was  investigated  which  in¬ 
creased  the  range  of  the  feedback  gain  directly.  The  specific  case  of  decreasing  the  detrimental 
effect  of  sensing  unmodelled  modes  is  thus  captured  as  a  special  case. 

Dynamic  Model  The  specific  pinned-pinned  beam  system  considered,  along  with  the  two 
sensor  domains,  is  shown  in  Figure  13.  A  single  point  force  actuator  was  postulated  as  the 
controlled  input,  located  a  distance  La  from  one  support. 
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The  dynamic  equation,  relating  external  loads  to  beam  vibration  is 


+  =  £(*)*«  (19) 

where  E  is  Young’s  modulus  and  /  is  the  area  moment  of  inertia  of  the  steel  beam,  p  is 
the  mass  density  and  A  is  its  cross-sectional  area,  and  L(x)  is  the  position  function  of  the 
actuator.  The  physical  dimensions  are  given  in  Table  1. 


Parameter 

Units 

Value 

La 

m 

0.06 

Lb 

m 

1.9 

Lsl 

m 

0.19 

LS2 

m 

0.38 

Wh 

m 

0.04 

tb 

m 

0.004 

P 

kg/m 3 

7800 

E 

N/m2 

210  xlO9 

p° 

e31 

m/V 

8.1987  xlO-2 

Table  1:  Physical  parameters  of  the  pinned-pinned  beam  of  the  example. 


An  assumed  modes  solution  takes  the  form 

n 

w(x,t)  =  ^2rn(t)4>i(x)  (20) 

i— 1 

where  r]i(t)  is  the  zth  generalized  coordinate  and  the  pinned-pinned  mode  shapes  (pi(x)  are 

4>i(x)  —  sin  (Pix)  (21) 

where  the  value  of  0i  can  be  determined  from  the  frequency  equation  of  the  beam,  sin(^Lb)  = 

0. 

In  general,  the  equation  of  motion  with  damping  can  be  written  as 

rji(t)  +  2CVA +  Arii{t )  =  Biu(t)  (22) 

where  £  is  the  matrix  with  modal  damping  coefficients,  A  is  the  matrix  of  eigenvalues  and 
the  control  weighting  coefficient  of  the  ith  mode,  Bi,  is 

Bi  =  J  <f>i{x)L(x)  dfl  (23) 

Cl 

where  fl  is  the  beam  domain.  The  state  space  representation  is 


x  =Ax  +  Bu 
y  =Cx 


(24) 
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(a)  Discrete  PVDF  sensors 


(b)  Segmented  distributed  PVDF  sensors 


Figure  13:  Side  and  top  views  of  the  pinned-pinned  beam  with  the  two  different  sensor  design 
domains  configurations  considered  in  the  design  example. 


where  the  2n  x  1  state  vector  x  contains  the  generalized  coordinates  and  their  speeds,  the 
2 n  x  2n  state  matrix  A  is 


A  = 


0 

-A 


I 

-2CVA 


(25) 


and  the  2n  x  1  input  weighting  vector  B  is 


(26) 


PVDF  Sensor  Equation  The  output  signal  was  assumed  to  be  measured  by  the  PVDF 
sensors,  shown  in  Figure  13.  The  model  of  the  charge  developed  on  the  fcth  piezoelectric 
sensor,  developed  by  Lee  and  Moon[13]  [14],  is 


Qk(t)  =  - 


Iff  d2w(x,y,t ) 

2  JZ{e31—dx^- 


+  e32 


d2w(x,y,t) 
dy 2 


) 


da:dy, 


where  z  =  (hp  +  hs)/2  with  hp  and  hs  being  the  thicknesses  of  the  piezoelectric  patch  and 
of  the  structure  respectively,  and  C31  and  632  are  piezoelectric  constants.  D*  represents  the 
domain  of  the  fcth  sensor,  while  w(x,  y,  t)  is  the  displacement  in  the  out-of-plane  direction. 
The  time  dependence  of  qk(t )  is  induced  because  of  the  time  dependence  of  w(x,y,t).  Sub¬ 
stituting  Eq.  (20)  into  Eq.  (27),  the  sensor  charge  is 


=  E 


i— 1 


ifz(e 31 


d2<j>i(x,y)  d24>j{x,y) 


dx 2 


+  632- 


dy 2 


0 


dxdy 


’h(t) 


(28) 
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Furthermore,  since  i(t)  —  dq/dt,  the  current  signal  due  to  the  mechanical  deformation  of  the 
structure  is 


1  f  (  d2<j>i(x,y)  d2<f>i(x,y)\ 

~2  JZ  {e31~~d^-  +  e32~W~ )  da^ 

£lk  . 


(29) 


In  practice,  if  displacement  is  the  measurement  desired,  a  charge  amplifier  should  be  con¬ 
nected  to  the  PVDF  patches,  and  if  the  measurement  desired  is  velocity,  a  current  amplifier 
should  be  used. 

The  optimization  process,  described  in  detail  in  the  next  section,  selects  the  sensor  dis¬ 
tribution  and  the  control  law  to  minimize  a  specified  objective  function.  To  accomplish  this, 
the  homogenization  approach  was  used  to  express  the  piezoelectric  coefficients,  e 3i  and  632, 
in  terms  of  the  optimization  parameters  av.  Specifically,  the  PVDF  material  was  modeled 
using  the  extended  unit  cell  microstructure  shown  in  Figure  14  which  was  described  in  de¬ 
tail  in  Section  5.1.  Each  unit  cell  is  characterized  by  the  square  void  with  edge  length  av. 
The  piezoelectric  coefficients  are  then  expressed  as  a  function  of  av  by  homogenizing  the 
cells  through  an  asymptotic  expansion,  taking  the  limit  to  an  infinite  number  of  holes  per 
unit  cell.  The  resulting  material  laws,  for  e3i(a)  and  e32(a„)  are  shown  in  Figure  15.  Note 
that  negative  values  of  av  are  permitted.  While  this  seems  inconsistent  from  the  physical 
description  of  a  ’’hole”,  mathematically  this  is  quite  convenient.  It  captures  the  notion  of 
oppositely  polled  PVDF  material  by  expressing  the  piezoelectric  coefficients  as 


4(a)  =  sign(a)4*(|  a  |)  (30) 

where  a  G  [—1,1],  and  e£*(|  a  |)  is  the  function  obtained  from  homogenization  [5]  valid  for 
o  G  [0, 1]. 


The  effect  of  the  parameterization  of  the  piezoelectric  coefficients  by  av  on  the  sensor 
equation,  Eq.  (28),  is  caputured  by  substituting  Eq.  (30)  into  Eq.  (28),  yielding  the  new 
sensed  charge  as 


»(*)  -  -\ /»{4,W*.»)]25^ + dxd y  (si) 

£4 
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Piezoelectric  stress  coefficient  e31  as  a  function  of  density  p 


Figure  15:  Homogeneous  as  a  function  of  density  p  for  a  microstructure 


The  output  matrix  C  of  Eq.  (24)  is  assembled  by  replacing  the  displacements  w(x,y,t )  by 
the  assumed  modes  expansion  of  Eq.  (20)  yielding  the  output  matrix  components 


cd(k,i)  =  - 

ilk 

.,d24>i{x, 

y)+4M*,v))d%£y) 

^  dxdy 

and  thus  the  output  matrix  is 

'Cd(  1,1) 

cd(  1,2)  . 

. .  0,(1,  n) 

0  0 

0 

C  = 

Cd(m,  1) 
0 

Cd(m,  2)  . 
0 

..  Cd(m,q) 
0 

0  0 

0,(1, 1)  0,(1, 2)  ... 

0 

0,(1,  n) 

0 

0 

0 

Cd(m,  1)  Cd(m,  2)  ... 

Cd(m,  q) 

(32) 


(33) 


For  the  one-dimensional  beam  example  considered  here  the  mode  shapes  <j>(x )  are  not  func¬ 
tions  of  y,  and  the  piezoelectric  coefficient  e\2  is  zero  thus  simplifying  the  calculation  of 
Eq.  (33). 


Control  Strategy  Description  In  the  system  simulation,  n  =  20  eigenmodes  were  con¬ 
sidered,  where  the  goal  was  to  damp  the  first  m  =  5  modes.  The  remaining  n  —  m  residual 
modes  were  left  in  the  model  to  capture  possible  spillover  effects.  The  control  design,  how¬ 
ever,  was  based  on  the  truncated  model  with  only  5  modes.  The  truncated  model  will  be 
denoted 


x  =Ax  +  Bu 
y  =Cx 


(34) 
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where  x  is  a  2m  x  1  truncated  state  vector,  A  is  a  2m  x  2m  state  matrix,  B  is  a  2m  x  1 
input  weighting  vector,  C  is  a  2 m  x  2m  output  matrix. 

A  state  feedback  controller  was  used  where 

u=  -Ktx  =  -KTRdx  (35) 

where  Rd  is  C~lC,  and  the  2m  x  1  gain  vector  K  was  obtained  using  the  LQR  method 

K  =  -BP  (36) 

r 

The  scalar  control  effort  weighting  is  denoted  r,  and  the  2m  x  2m  matrix  P  is  the  solution 
to  the  matrix  Ricatti  equation 

AtP  +  PA--PBBtP  +  Q  =  0  (37) 

r 

where  the  weighting  matrix  for  the  state  Q  was  chosen  as 


Q  = 


A  0 
0  I 


(38) 


Based  on  r  and  Q,  the  average  energy  (both  kinetic  and  potential)  contained  in  the  oscillation 
and  total  energy  consumed  by  the  controller  are  minimized  at  the  same  time.  This  minimum 
value  of  r  reflects  the  maximum  feedback  gain  resulting  in  a  neutrally  stable  system. 

Although  LQR  is  used  in  this  study,  it  should  be  noted  that  in  practice  the  required  state 
vector  is  not  measured  directly,  and  thus  a  state  estimator  would  need  to' be  used  [20].  This, 
however,  does  not  change  the  main  results  of  this  work  and  the  focus  on  combined  sensor 
and  control  law  optimization. 


Formulation  of  the  Optimal  Design  Problem  Simultaneous  controller  and  system 
optimization  problems  can  be  classified  into  one  of  two  types,  (1)  Direct  Gain  Optimization 
or  (2)  Indirect  Gain  Optimization.  The  direct  gain  optimization  approach  assumes  a  postu¬ 
lated  control  strategy  with  some  gains  or  parameters  unspecified.  A  numerical  optimization 
process  is  then  used  to  select  the  gains  of  the  controller  and  the  parameters  that  define  the 
structure’s  design  to  minimize  some  performance  criteria.  The  indirect  approach,  as  em¬ 
ployed  in  this  work,  seeks  to  enhance  the  system’s  performance  without  explicitly  specifying 
the  actual  gains  of  the  controller.  In  this  example,  reducing  spillover  sensitivity  is  the  main 
goal,  for  a  specified  controller  form  (LQR).  Once  the  sensor  topology  is  formed  through  the 
optimization  process,  the  designer  has  the  ability  to  select  a  wide  range  of  possible  controller 
gains.  This  method  has  the  attractive  feature  of  yielding  a  controller /system  design  that  is 
more  flexible  as  compared  to  the  direct  gain  approach.  In  the  remainder  of  this  section  the 
optimization  problem  is  formulated  for  two  different  cost  functions. 

Each  PVDF  sensor,  for  both  the  discrete  and  the  segemented  distributed  design  domains 
of  Figure  13,  was  divided  into  p  segments.  The  piezoelectric  properties  of  each  segment  was 
uniquely  characterized  by  its  density-like  parameter  av.  The  set  of  5 p  parameters  (since  there 
are  5  sensors),  denoted  A,  constituted  the  only  free  variables  in  the  optimization  process. 
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Two  cost  functions  were  considered.  The  first  one,  denoted  J\,  was  based  on  the  observ¬ 
ability  gramian,  where  the  idea  was  to  improve  the  system’s  stability  margin  by  suppressing 
the  observation  spillover  of  the  residual  modes.  The  other,  denoted  J 2,  was  simply  the  LQR 
control  weighting  scalar  r. 

According  to  Hac  (1993)[11),  if  the  natural  frequencies  of  the  structure  are  well  spaced 
and  the  damping  coefficients  are  small,  then  the  system’s  observability  gramian  can  be 
approximated  as 


W0  = 


4C2^2 


0 


(39) 


The  values  Q  and  w*  are  the  modal  damping  coefficients  and  natural  frequencies  of  the  open 
loop  system.  The  quantity  Pa  is  defined  as 


Pa  =1  if  i  €  l..m 

m 

Pu='Y^Rl{i,k)  if  iem  +  l..q 
fc= 1 

Using  this  notation,  the  first  cost  function  is  defined  as 


(40) 


where  V  is  the  set  of  modes  whose  observability  gramian  entries  are  to  be  minimized.  Since 
the  goal  was  to  control  the  first  m  modes,  V  should  contain  the  modelled  modes  greater 
than  m.  The  superscript  c  on  and  Q  indicates  a  closed-loop  quantity.  A^fn  and  A^  are 
minimum  and  maximum  eigenvalues  of  Cd ■  The  term 


denotes  the  damped  conditioning  number.  This  definition  is  computationally  meaningful 
only  provided  the  conditioning  number  of  the  matrix  (the  ratio  of  the  largest  to  the  smallest 
eigenvalue)  is  not  larger  than  the  inverse  of  the  machine  accuracy  or  relative  round  off  error. 
For  the  smart  structure  application,  we  desire  good  sensing  of  the  first  m  modes,  which  the 
conditioning  number  of  Cd  ensures.  77  denotes  a  numerical  damping  coefficient  [23]  that  is 
selected  in  an  ad-hoc  manner  to  aid  the  optimization  convergence.  In  this  study  a  value  of 
77  —  1.6  was  used. 

The  optimization  problem  can  be  summarized  as  finding  the  parameter  set  A  that  mini¬ 
mizes  Eq.  (41)  whose  terms  Pa,  (f,  and  u>?  are  the  result  of  implementing  an  LQR  controller, 
designed  from  the  m  mode  truncated  model  of  Eq.  (34),  on  the  n  mode  design  model  of 
Eq.  (24)  .  Each  design  parameter  of  A  is  constrained  such  that  0  <  a*  <  1. 
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The  second  cost  function  is  the  LQR  control  weighting  term,  or, 

J2(a)  =  min(r)  (43) 

with  the  constraint  that  all  closed-loop  poles  have  the  property  that 

Reijpi)  <0  (44) 

This  ensures  the  system  is  stable.  Small  values  of  the  cost  function  J2{a)  mean  that  large 
values  of  the  overall  closed-loop  gain  can  be  permitted  without  causing  instability. 

Similar  to  the  previous  cost  function,  the  optimization  problem  can  be  summarized  as 
finding  the  parameter  set  A  that  minimizes  Eq.  (43)  such  that  the  closed-loop  poles  of  the 
n  mode  design  model  of  Eq.  (24)  are  stable  where  its  LQR  controller  is  designed  based  on 
the  m  mode  truncated  model  of  Eq.  (34).  Again,  each  design  parameter  of  A  is  constrained 
such  that  0  <  a*  <  1. 

A  sequential  quadratic  programming,  numerical  optimization  code  was  used  to  solve  both 
optimization  problems  where  the  initial  values  of  all  the  design  parameters  were  set  to  1.0. 
Due  to  the  complexity  of  the  problem,  no  proof  was  found  to  indicate  a  globally  optimal 
solution.  From  a  design  perspective,  this  is  acceptable  since  the  optimization  process  does 
provide  a  better  design  than  using  the  nominal,  full-density  sensor  distribution.  Finally,  since 
the  pinned-pinned  beam  mode  shapes  are  symmetric  about  the  geometric  center,  A  half-beam 
optimization  approach  was  used,  thus  reducing  the  number  of  optimizeable  parameters  by 
one  half. 


Initial  Distribution  Jx  Distribution  J2  Distribution 
~JX  121.68  0.9368  116.8091 

J2  0.0594  0.0103  0.00385 _ 


Table  2:  Cost  Function  Values  of  Optimal  Distribution  of  Discrete  PVDF 


Initial  Distribution  J*  Distribution  J2  Distribution 
~JX  3.0524  0.0468  0.6887 

J2  0.0267  0.00073  0.000249 


Table  3:  Cost  Function  Values  of  Optimal  Distribution  of  Segmented  Distributed  PVDF 


Optimal  Sensor  and  Controller  Design  Results  For  each  sensor  configuration  the 
sensor  density  and  control  law  were  optimized  simultaneously  for  the  two  cost  functions 
described  in  the  previous  section.  Table  2  summarizes  the  results  for  the  discrete  sensor 
case,  and  Table  3  shows  the  results  for  the  segmented,  distributed  case.  In  addition  to 
showing  the  change  in  cost  from  the  initial  densities  (pv  =  1)  to  the  optimal,  the  non 
optimized  cost  function  is  evaluated.  For  example,  the  first  row  of  Table  2  shows  the  initial 
value  of  Ji,  the  value  of  Jx  (where  *  denotes  an  optimal  quantity)  and  finally  the  value  of 
J\  using  the  distribution  that  produced  J72*.  This  helps  to  illustrate  the  differences  between 
cost  functions. 
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0  0.2  0.4  0.6  0.6  1  1.2  1.4  1.6  1.8  2 


x(m) 


(a)  Beam  with  the  discrete  PVDF  sensor  configuration 
(,7i=0.9368) 


0  0.2  0.4  0.6  0.8  1  1.2  1.4  1.6  1.8  2 


x(m) 

(b)  Beam  with  segmented,  distributed  PVDF  sensor 
configuration  (Ji  =0.468) 

Figure  16:  The  optimal  normalized  PVDF  density  as  a  function  of  beam  position  for  both  sensor 
configurations.  In  both  cases  the  PVDF  was  optimized  using  the  observability  gramian 
cost  function 
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(a)  Beam  with  the  discrete  PVDF  sensor  configuration 
(*72=0.00385) 


0  0.2  0.4  0.6  0.8  1  1.2  1.4  1.6  1.8  2 


x(m) 

(b)  Beam  with  segmented,  distributed  PVDF  sensor 
configuration  (72=0. 000 249) 

Figure  17:  The  optimal  normalized  PVDF  density  as  a  function  of  beam  position  for  both  sensor 
configurations.  In  both  cases  the  PVDF  was  optimized  using  the  LQR  stability  margin 
cost  function  7b- 
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(a)  Beam  with  the  discrete  PVDF  sensor  configura¬ 
tion 
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(b)  Beam  with  the  segmented,  distributed  PVDF  sen¬ 
sor  configuration 

Figure  18:  Pole  locations  with  uniform  density  distributions.  Mode  numbers  are  given  for  select 
mode  shapes. 
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a)  Beam  with  the  approximate  discrete  PVDF  sensor 
onfiguration  (^2=0.0046) 


x  (m) 


(b)  Beam  with  the  approximate  segmented,  distributed 
PVDF  sensor  configuration  (J2 =0.000965) 


Figure  19:  The  approximate  optimal  normalized  PVDF  density  for  both  sensor  configurations. 
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Observability  Gramian  Values 

Mode  Initial  J\  Optimal  J2  Optimal 


Number 

Disc. 

Seg.  Dist. 

Disc. 

Seg.  Dist. 

Disc. 

Seg.  Dist. 

6 

0 

0.1908 

0.0129 

0.0464 

0.0949 

0.1275 

7 

0.1939 

0.3392 

0.2862 

0.0457 

1.8330 

0.3320 

8 

0.4995 

0.7631 

0.0033 

0.0010 

0.0053 

0.0027 

9 

1.6053 

3.0524 

0.8714 

0.0350 

1.5881 

0.5292 

10 

7.9914 

0 

0.9160 

0.0275 

9.8149 

0.4927 

11 

121.6804 

3.0524 

0.9368 

0.0430 

116.8091 

0.6887 

12 

0 

0.7631 

0.9614 

0.0468 

1.5047 

0.3930 

13 

99.0245 

0.3392 

0.9614 

0.0466 

8.1072 

0.2689 

14 

5.2304 

0.1908 

4.2907 

0.4026 

1.0705 

0.2266 

15 

0.8228 

0.1221 

22.7815 

1.5989 

14.7049 

0.3220 

16 

0.1908 

0.1908 

4.0236 

5.7998 

2.0425 

0.4112 

17 

0.0503 

0.3392 

3.4876 

50.2938 

2.2847 

1.5027 

18 

0 

0.7631 

0.1794 

75.7091 

1.5628 

1.1079 

19 

0.006 

3.0524 

1.8262 

922.4670 

11.1493 

1.5256 

20 

0 

0 

7.1655 

1.1942 

0.3801 

0.0331 

Table  4:  Comparison  of  observability  gramian  values  for  both  sensor  domains  (‘Disc*  =  Discrete, 
‘Seg.  Dist’  =  Segmented  Distributed).  In  addition  to  the  nominal  distributions,  both 
the  J\  and  J2  optimal  distributions  are  shown. 


Discussion  When  considering  the  lower  values  of  Table  3  it  appears  that  regardless  of 
the  cost  function  used,  the  segmented,  distributed  sensor  configuration  yields  better  results 
than  that  of  the  discrete  sensor.  The  only  caveat  to  this  conclusion  is  illustrated  in  the 
4th  column  of  Table  4  which  shows  the  observability  gramians  for  each  uncontrolled  mode 
of  the  simulation.  Specifically,  when  the  segmented,  distributed  sensor  is  used  with  the 
observability  gramian  cost  function,  J\ ,  residual  modes  17-19  acquire  a  large,  and  potentially 
poor,  observability  gramian  value.  It  should  be  noted  that  in  the  last  column  of  Table  4  does 
not  show  the  opposite  effect.  That  is,  when  the  distribution  is  optimized  with  respect  to 
the  LQR  stability  margin  the  observability  gramian  values  do  not  have  any  extreme  values. 
These  features  suggest  that  the  segmented  distributed  sensor  configuration,  along  with  LQR 
stability  margin  optimization,  J2,  is  the  best  choice. 

This  conclusion  is  further  justified  when  considering  the  smoothness  of  the  actual  PVDF 
density  profiles  of  Figure  16  and  Figure  17.  Specifically,  fabrication  of  variable  density 
active  material  is  a  rather  new  discipline  with  limitations  on  gradients.  In  all  cases  the  LQR 
stability  margin  optimal  distributions  were  smoother  than  the  observability  gramian  ones. 
This  is  especially  true  for  the  segmented  distributed  configuration  of  Figure  17(b)  where  the 
density  profile  may  be  approximated  by  a  simple  trapezoidal  shape. 

It  can  be  seen  in  Figure  16  and  Figure  17  that  sensor  density  distribution  depends  on 
what  cost  function  is  used.  The  rmin  optimal  distribution  of  Figure  17  has  all  positive 
densities.  This  can  be  interpreted  as  all  sensors  being  on  one  side  of  the  beam.  This  is 
an  important  result,  because  it  implies  that  better  performance  of  the  structure  is  achived 
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when  sensors  axe  distributed  uni-sided.  Practically,  it  could  imply  simpler  manufacturing 
and  more  cost-effective  active  structures. 

The  shape  of  the  distribution  also  has  some  physical  significance.  To  illustrate  this,  we 
will  focus  on  the  rmin  optimal  cases  of  Figure  17.  The  goal  of  the  rmin  optimization  is  to  find 
the  parameter  set  A  that  allows  for  the  largest  LQR  gains  while  still  maintaining  stability. 
Increasing  the  gains  beyond  this  point  results  in  closed-loop  poles  migrating  to  the  right  half 
of  the  complex  plane.  The  shape  of  the  optimal  distribution  is  related  to  the  mode  shapes 
of  the  poles  that  are  the  first  to  go  unstable.  To  illustrate  this,  plots  were  constructed  that 
show  the  open  loop  and  closed-loop  poles  for  both  the  discrete  and  segmented  distributed 
sensor  configurations,  shown  in  Figure  18.  For  the  discrete  sensor  case  of  Figure  18(a), 
modes  7,  13,  14,  and  15  are  most  likely  to  go  unstable  as  the  loop  gain  is  increased.  For 
the  segmented  distributed  case  of  Figure  18(b),  modes  8,  9,  11  and  12  go  unstable.  The 
density  distributions  of  Figure  17(a)  and  Figure  17(b)  are  strongly  influenced  by  these  mode 
shapes.  More  specifically,  approximate  optimal  solution  can  be  constructed,  and  are  shown 
in  Figure  19,  where  the  distributions  are  formed  as  linear  combinations  of  these  mode  shapes. 
The  equations  used  to  generate  the  approximate  densities  of  Figure  19  are 


Pa  =  3  -I-  07  —  0.3013  —  014  +  015  (45) 

Pb  =  10  —  0.508  —  09  —  011  —  0.3012  (46) 

where  the  mode  shape  0*  are  denoted  in  Eq.  21,  and  the  suboptimal  rmin  cost  functions  are 
0.0046  and  0.000965  for  the  densities  of  Figure  17(a)  and  Figure  17(b)  respectively.  This 
compares  favorably  to  the  optimal  costs  of  0.00385  and  0.000249.  To  summarize,  the  shape 
of  the  optimal  distributions  allow  the  controller  to  better  sense  those  modes  that  contribute 
the  most  to  limiting  the  stability  margin. 

5.3.2  Simultaneous  Optimization  of  an  Active  Structure 

The  main  objective  of  this  research  is  to  develop  and  assess  a  method  for  simultaneous  opti¬ 
mization  of  active  structures  considering  the  conventional  structure  topology,  active  material 
layout,  and  the  control  system.  In  this  section  this  method  is  used  to  demonstrate  the  ro¬ 
bust  control  of  a  piezoelectric  laminate  beam  by  simultaneously  optimizing  the  conventional 
material  distribution,  active  material  distribution  and  the  control  system  for  a  closed-loop 
system.  The  parallel  genetic  algorithm  based  active  structure  design  code  was  used,  includ¬ 
ing  ABAQUS  for  model  development  and  Matlab®  for  control  design.  The  structure  was 
initialized  as  a  grid  of  cells  with  complete  structural  and  energetic  properties.  During  the 
optimization  process,  structural  material  and  active  material  was  removed  from  the  cells 
to  minimize  a  multi-objective  cost  function.  In  order  to  demonstrate  the  improvement  of 
the  control  efficiency,  the  optimally  designed  system  was  compared  to  a  conventional  design 
consisting  of  a  cantilever  alluminum  beam  with  a  PZT  patch  actuator  at  the  root  and  2 
PVDF  sensors  fully  covering  the  beam’s  underside. 

Dynamics  of  a  Piezoelectric  Laminate  Beam  The  cantilever  flexible  structure  with 
one  PZT  actuator  and  two  PVDF  sensors  as  shown  in  Figure  20  was  considered.  The  physical 
dimensions  for  the  beam  are  given  in  Table  5. 
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Parameter 

Units 

Value 

beam 

U 

m 

0.32 

Wb 

m 

0.03 

tb 

m 

0.001 

p 

kg/m3 

2700 

E 

N/m2 

69  xlO9 

PZT 

La 

m 

0.8 

wa 

m 

0.03 

ta 

m 

0.267  xlO-3 

Xi 

m 

0 

X2 

m 

0.08 

e31  a 

m/V 

12.54 

e32  a 

m/V 

12.54 

PVDF 

Ls 

m 

0.32 

Ws 

m 

0.019 

ts 

m 

9  xlO-6 

P° 

e31s 

m/V 

0.046 

P° 

e32s 

m/V 

0.046 

Table  5:  Physical  parameters  of  the  cantilever  beam  of  the  example. 


u 


Figure  20:  Cantilever  beam  model. 
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The  beam’s  transverse  deflection  at  point  x  and  time  t  is  denoted  w(x,  t),  assuming 
the  beam  as  a  one-dimensional  system  only.  The  partial  differential  equation  (PDE)  which 
describes  the  dynamics  of  the  beam  is 

r,Td4w(x,t)  Ad2w(x,t)  d2Ma(x,t)  . 

EI~W~  +  pA~dP~  =  W  ~  (4?) 

where  p,  A,  E,  and  I  represent  density,  cross-sectional  area,  Young’s  modulus  of  elasticity, 
and  moment  of  inertia  about  the  neutral  axis  of  the  beam,  respectively. 

The  general  equation  of  motion  can  be  described  as 

Vi(t)  +  2C^iji(t)  +  Ar)i(t)  =  BiVa(t)  (48) 


where  C  are  the  modeal  damping  coefficients,  A  is  the  matrix  of  eigenvalues  and  the  control 
weighting  coefficient  of  the  ith  mode,  Bi,  is 


cP&jxpy) 
dx 2 


+  ^32  a 


d24>j(x,y) 

dy 2 


) 


dxdy 


(49) 


where  fl  is  the  actuator  domain.  Using  the  homogenization  approach  to  map  cell  densities  to 
material  properties  and  Eq.  (28),  the  sensor  equation  and  the  actuator  equation  in  response 
to  density  can  be  obtained 


Cd(k,i) 


Bi  = 
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(50) 
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The  objective  of  the  experiment  was  to  increase  the  stability  margin  of  the  closed-loop 
system  and  control  the  beam’s  first  two  modes  of  vibration  using  an  LQR  control  law  with 
the  cost  function  of  Eq.  43.  A  block  diagram  of  the  closed-loop  system  is  shown  in  Figure  21 
Densities  of  the  PVDF  sensor,  the  PZT  actuator  and  the  beam  were  considered  optimizeable 


Figure  21:  Control  block  diagram. 

parameters  of  the  cost  function  Rmin  which  provides  a  wide  range  of  possible  gains  to  ensure 
the  system  is  stable.  Simply,  with  optimal  density  distributions,  the  designer  has  the  ability 
to  select  a  larger  gain  in  comparison  to  that  of  a  nominal,  uniform  density  design. 
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Simulation  and  Experimental  Results  The  plant  model  used  during  optimization  cap¬ 
tures  the  first  3  modes  of  the  beam.  Since  the  control  objective  focuses  on  only  the  first 
two  modes,  the  third  mode  is  included  so  that  the  final  design  is  more  resiliant  to  third 
mode  spillover.  To  illustrate  the  spillover  effect  for  the  nominal,  uniform  density  beam,  its 
closed-loop  pole  locations  are  shown  in  Figure  22.  From  this  figure  it  is  evident  that  the 
system  will  eventually  become  unstable  as  the  feedback  gain  increases. 

Before  beginning  the  optimization  process  it  is  important  that  the  model  be  predictive 
of  the  physical  system.  To  gauge  this,  the  nominal  beam  with  PVDF  sensors  was  tested  and 
compared  to  a  simulation  of  the  3  mode  model.  The  resulting  comparison  of  this  open-loop 
case  is  shown  in  Figure  23. 
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Figure  22:  Pole  locations  of  the  system  with  initial  density  setup. 

The  objective  of  the  optimization  process  was  to  find  a  design  defined  by  the  density 
distributions  of  the  conventional  beam,  PZT  actuator,  and  PVDF  sensors  yielding  a  larger 
stability  margin  than  the  nominal  uniform  density  case.  The  design  environment  described 
earlier,  consisting  of  a  parallel  processing  Genetic  Algorithm  opimization  engine,  ABAQUS 
for  model  generation,  and  Matlab®  for  control  design,  was  used  with  20000  members  in  the 
initial  population.  20  population  generations  were  created  where  the  simultaneous  design 
code  ran  on  20  workstations  in  parallel.  The  ’’raw”  optimal  solution  is  shown  in  Figure  24(a), 
and  yielded  an  rmin  of  0.004  as  compared  to  rmin  =  0.3165  for  the  nominal  beam  with 
completely  covered  PVDF  and  PZT  material.  Since  the  discontinuous  distributions  would 
be  difficult  to  machine,  smoothed  distributions  were  generated  using  a  6th  order  polynomial 
fit  to  the  optimal  distributions,  similar  to  the  approach  used  in  Section  5.3.1.  This  changed 
the  cost  only  slightly,  increasing  rmin  to  0.0041.  These  are  plotted  in  Figure  24(b).  It 
should  be  noted  that  the  widths  of  the  various  materials  were  constrained  according  to: 
PZT  material  >  0.5,  PVDF  material  >  0.5  and  the  conventional  material  >  0.65. 

Figure  25  shows  time  histories  of  the  nominal  closed-loop  system  to  the  optimal  closed- 
loop  system.  Since  the  rmin  was  decreased  by  a  factor  of  80  for  the  optimal  system,  the  loop 
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(a)  Signals  of  the  1st  PVDF  sensor 


(b)  Signals  of  the  2nd  PVDF  sensor 

FIGURE  23:  Simulation  and  Experiment  sensor  signals  of  the  open  loop  system. 


5  ACCOMPLISHMENTS 


45 


0  0.05  0.1  0.15  0.2  0.25  0.3  0.35 

length  of  sensor  (m) 


0  0.05  0.1  0.15  0.2  0.25  0.3  0.35 

length  of  beam  (m) 


0  0.01  0.02  0.03  0.04  0.05  0.06  0.07  0.08  0.09  0.1 

length  of  actuator  (m) 


(a)  Optimal  distributions  using  GA  method 
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(b)  Smooth  optimal  distributions  after  interpolated 
Figure  24:  Optimal  distributions  for  conventional  beam,  PZT  actuator  and  PVDF  sensor. 
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gain  could  be  increased  as  compared  to  the  nominal  system  while  still  ensuring  closed-loop 
stability.  This  resulted  in  a  minimum  decrease  in  strain  by  a  factor  of  1/2  as  can  be  seen  in 
Figure  25 


(a)  Signals  of  the  1st  PVDF  sensor 


(b)  Signals  of  the  2nd  PVDF  sensor 

Figure  25:  Simulated  sensor  signals  of  the  optimal,  closed-loop  system. 
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6  Publications 

The  publications  listed  below  are  divided  into  those  in  archival  journals,  conference  proceed¬ 
ings,  and  theses. 

6.1  Peer-Reviewed  Publications 

1.  Buehler,  M.  J.,  Bettig,  B.  P.,  and  Parker,  G.  G.,  2003,  ‘Numerical  Homogenization 
of  Smart  Material  Finite  Element  Cells,’  Communications  on  Numerical  Methods  in 
Engineering ,  Accepted  for  Publication. 

2.  Chen,  W.,  Buehler,  M.,  Parker,  G.  G.,  and  Bettig,  B.,  2003,  ‘Optimal  Sensor  Design 
and  Control  of  Piezoelectric  Laminate  Beams,’  IEEE  Transactions  on  Control  System 
Technology,  Accepted  for  Publication. 

3.  Buehler,  M.  J.,  Bettig,  B.  P.,  and  Parker,  G.  G.,  2003,  ‘Topology  Optimization  of 
Smart  Structures  using  a  Homogenization  Approach’,  Journal  of  Intelligent  Material 
Systems  and  Structures,  Accepted  for  Publication. 

6.2  Publications  in  Conference  Proceedings 

1.  Chen,  W.,  Buehler,  M.,  Parker,  G.  G.,  and  Bettig,  B.,  2003,  ‘Robust  Design  and 
Control  of  Piezoelectric  Laminate  Beams  Using  a  Simultaneous  Optimization  Method,’ 
SPIE  10th  Annual  International  Symposium  on  Smart  Structures  and  Materials,  March 
3-6. 

2.  Buehler,  M.,  Bettig,  B.,  and  Parker,  G.  G.,  2002,  ‘Topological  Optimization  of  Smart 
Structures  Using  an  Homogenization  Approach,’  SPIE  9th  Annual  International  Sym¬ 
posium  on  Smart  Structures  and  Materials,  San  Diego  CA,  March  18-20. 


6.3  Theses 

1.  Buehler,  M.,  2001,  ‘Homogenization  of  Smart  Material  Cells  for  Topological  Optimiza¬ 
tion,’  Master’s  Thesis,  Michigan  Technological  University. 
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7  Interactions  and  Transitions 

Professors  Parker  and  Bettig  have  communicated  the  capabilities  of  the  Parallel  Genetic 
Algorithm  Active  Structure  Design  Code  to  contacts  at  AFRL/MNAV.  We  are  currently- 
exploring  their  interest  in  using  it  for  cruise  missile  active  wing  twist  design.  This  would 
require  some  new  development  in  the  area  of  implementing  aerodynamic  load  models  into 
the  existing  finite  element  portion  of  the  code. 

Another  notable  interaction  has  developed  through  one  of  the  M.S.  students  working  on 
the  project.  Specifically,  Markus  Buehler  received  his  M.S.  degree  from  Michigan  Techno¬ 
logical  University  in  2001,  and  then  went  to  the  Max-Planck  Institute  in  Germany.  He  is 
still  an  active  participant  in  the  project,  and  has  facilitated  important  technical  relationships 
between  Michigan  Tech  and  the  Max-Planck  Institute. 
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8  New  Discoveries,  Inventions,  and  Patent  Disclosures 

The  project  resulted  in  the  creation  of  the  Parallel  Genetic  Algorithm  Active  Structure 
Design  Code  described  in  detail  in  the  Accomplishment  section  of  the  report.  Since  it 
relies  on  a  general,  off-the-shelf  structural  analysis  code  (ABAQUS)  and  an  equally  powerful 
control  design  tool  (Matlab®  ),  its  applicability  is  broad  and  covers  most  active  structure 
design  possibilities. 
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9  Honors  and  Awards 

•  G.  Parker  was  promoted  to  the  rank  of  Associate  Professor  in  2001. 

•  G.  Parker  received  the  Michigan  Technological  University  Distinguished  Teaching  Award 
in  2001. 

•  G.  Parker  received  the  Society  of  Automotive  Engineers  Ralph  E.  Teetor  Award  in 

2002. 
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